Astronomy &: Astrophysics manuscript no. 3439 


February 5, 2008 


(DOI: will be inserted by hand later) 





Probing Turbulence with Infrared Observations in OMCl 

M.Gustafsson^*, D.Field^, J.L.Lemaire^, and F.P.Pijpers^ 

^ Department of Physics and Astronomy, University of Aarhus, DK-8000 Aarlius C, Denmark 

e-mail: maikeng@phys .au.dk 
^ Observatoire de Paris & Universite de Cergy-Pontoise, LERMA & UMR 8112 du CNRS, 92195 Meudon, France 
^ Space and Atmosplieric Pliysics, Dept. Physics, Imperial College London, England 

For Main Journal, Diffuse Matter in Space. Received: Accepted: 

Abstract. A statistical analysis is presented of the turbulent velocity structure in the Orion Molecular Cloud 
at scales ranging from 70 AU to 3-10* AU. Resuhs are based on IR Fabry-Perot interferometric observations of 
shock and photon-excited H2 in the K-band S(l) v=l-0 line at 2.121/im and refer to the dynamical characteristics 
of warm perturbed gas. Data consist of a spatially resolved image with a measured velocity for each resolution 
limited region (70AUx70AU) in the image. The effect of removal of apparent large scale velocity gradients is 
discussed and the conclusion drawn that these apparent gradients represent part of the turbulent cascade and 
should remain within the data. Using our full data set, observations establish that the Larson size-linewidth 
relation is obeyed to the smallest scales studied here extending the range of validity of this relationship by nearly 
2 orders of magnitude. The velocity probability distribution function (PDF) is constructed showing extended 
exponential wings, providing evidence of intermittency, further supported by the skewness (third moment) and 
kurtosis (fourth moment) of the velocity distribution. Variance and kurtosis of the PDF of velocity differences 
are constructed as a function of lag. The variance shows an approximate power law dependence on lag, with 
exponent significantly lower than the Kolmogorov value, and with deviations below 2000AU which are attributed 
to outflows and possibly disk structures associated with low mass star formation within OMCl. The kurtosis 
shows strong deviation from a gaussian velocity field, providing evidence of velocity correlations at small lags. 
Results agree accurately with semi-empirical simulations in Eggers & Wane ( 199^. 

In addition, 170 individual H2 emitting clumps have been analysed with sizes between 500 and 2200 AU. These 
show considerable diversity with regard to PDFs and variance functions (related to second order structure func- 
tions) displaying a variety of shapes of the PDF and different values of the scaling exponent within a restricted 
spatial region. However, a region associated with an outfiow from a deeply embedded O-star shows high values of 
the scaling exponent of the variance function, representing a strong segregation of high and low exponent clumps. 
Our analysis constitutes the first characterization of the turbulent velocity field at the scale of star formation and 
provide a dataset which models of star-forming regions should aim to reproduce. 

Key words. ISM: individual objects: OMCl, ISM: kinematics and dynamics, ISM: molecules, shock waves, turbu- 
lence, infrared:ISM 



1. Introduction 

A key to understanding the mechanism of star forma- 
tion is to characterise in detail the nature of the turbu- 
lent, weakly ionized and magnetised plasma in which stars 
form. Recently Gustafsson et al. 2003 (Paper I) published 
observational results for the vibrational emission of H2 in 
the archetypal star-forming region in the Orion Molecular 
Cloud, OMCl. Using a combination of Fabry-Perot intcr- 



* Based on observations obtained at the Canada-France- 
Hawaii Telescope (CFHT) which is operated by the National 
Research Council of Canada, the Institut National des Science 
de rUnivers of the Centre National de la Recherche Scientifique 
of France, and the University of Hawaii. 



ferometry and the PUEO adaptive optics system on the 
Canada- France-Hawaii Telescope (CFHT), data achieved 
a spatial resolution of O^.'IS (70AU at the distance of 
OMCl, D=460 pc dBallv et all bOOOl) ). with a velocity 
discrimination of 1 kms~^ in regions of high brightness. 
These data, limited as they are to highly excited regions 
in OMCl, provide the first opportunity to study the phys- 
ical properties of a star-forming region both in terms of 
its morphology and bulk gas motion at the scale of star 
formation. The aim of the present paper is to characterise 
the nature of the turbulent velocity field in OMCl for the 
subset of regions represented through vibrationally excited 
H2. Our results should help to provide a benchmark for 
comparison with MHD models of star-forming regions. 
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The importance of turbulent gas motion in regulating 
star formation may be illustrated as follows. Relative bulk 
motions obtained from observations reported here range 
up to 40 kms~^ in dense gas. Since these motions are su- 
personic, it is evident that bulk motion must contribute 
more energy per unit volume, that is, pressure, than ther- 
mal energy. Equally if a simple scaling law is used between 
magnetic induction, B, and par ticle density, n, whereb y 
B=bni/VG, with n in cm'^ l|Trola,nd Heilej ll 98(il) . 
then the energy contained in bulk mass motion will ex- 
ceed the energy in the magnetic field, (in SI units), 
for bulk velocities greater than 2. 2 b kms~^. The constan t 
b lies typically between 1 and 2 l)Kristensen et alJl2005|) . 
Thus in many regions the turbulent pressure may exceed 
the magnetic pressure, given that bulk motion proceeds at 
tens of kms~^as mentioned above. The magnetic pressure 
in turn typically exceeds the thermal pressure, the ratio of 
thermal to magnetic pressure being ^5x10""^ T/b^, where 
T is the gas temperature. In the above simple prescription, 
any anisotropy of the magnetic field has been ignored. At 
all events, on purely energetic grounds, in earlier stages of 
star formation, turbulence is the controlling support mech- 
anism against gravitational collapse, if bulk mass motions 
exceed a few kms~^. This statement remains true even in 
the warm regions involved in the present study. 

Observational and computational evidence for the sim- 
ple picture presented above has grown consider ably in 
the la st de cade as described in recent reviews bv iLarsonl 
l)2003(l and lMac Low fc Klessei] l|2004 . Moreover the role 
of turbulence goes further than simply slowing the pro- 
cess of gravitational collapse. Simulations show that tur- 
bulence may determine the initial mass function (IMF) 
for star formation through turbulent fragmentation (e.g. 
iNordlund fc Padoan 20031 It is evident that MHD mod- 
els, which make such far-reaching predictions, need to be 
constrained by observations. The meeting of theory and 
observation may be achieved by recording the statistical 
properties of velocity fields and this standard approach 
is adopted here. In the present work we make only very 
limited comparison with published models. Rather, we 
present modellers with characteristic behaviour which it 
is their challenge to reproduce. 

Several techniques have previously been used to char- 
acterize the structure of brightness and velocity in molec- 
ul ar clouds based largely upon CO observat i ons, r eviewed 
bv lGoodman et al. and lMiesch et alJ l|l999l) . These 

earlier observations, tracing relatively cool and low density 
gas, are limited in spatial resolution and can only be used 
to address the physics at scales larger than roughly O.OSpc 
(6000 AU). When dealing with turbulence it is assumed 
that most features are scale-free between the driving and 
dissipation scale and that the behaviour for example of 
the structure functions (see Sect. l4.4|l can be extrapolated 
to the very much smaller scales of individual star forma- 
tion. This is however far from certain. The high spatial 
resolution infrared data reported in Paper I give the op- 
portunity to discover how structure functions develop at 
much smaller scales than have hitherto been probed. 



In this paper we use three standard statistical mea- 
sures to provide observational constraints on theories 
and simulations concerning gravitational collapse and star 
formation in a turbulent ISM. First, we test whether 
the observed scaling behaviour of the velocity structure, 
encaps ulated i n the size-linewidth relation or "Larson 
law" l|Larsonl Il98ltl . holds for the smaller scales in- 
volved here (Sect. 14. 2|) . Second, we use probabihty dis- 
tribution functions (PD Fs) of velocities (Sec t. I4.3[l as a 
probe of interrnittency j Palgarone fc Phillipslil99 flL .19911: 
Fakarone et al.lll994l: iMiesch fc Scalolll995llMiesch et alJ 



19991 lOssenkopf fc Mac Lowl l2002l) . Thhd, we calcu- 
late t he low order moments of the velocity difference 
PDF (Miesch fc Sca,lfjri99ilLis et a,]jri99a IMiesch ^TZI 
'1999; 'Os senkonf fc Mac Lowlliool as a function of lag 
(Sect. I4.4|l . that is, the structure functions or functions 
directly related to structure functions. In addition, the 
high spatial resolution of the data allows us to go further 
than determining these quantities for the observed field as 
a whole. In Sect. |H|we identify 170 individual clumps in 
OMCl with a mean size of 1300 AU. We calculate the indi- 
vidual PDFs of peak velocities and the variance functions 
for these clumps. 

2. Observations and data reduction 

Near infrared K-band observations of the strongest H2 
emission line in Orion, v= 1-0 S(l) at 2.121 /im, are used 
as a tracer of radial gas velocity to provide the data for 
the statistical analysis performed here. The observational 
data are the same as those described in Paper I and a 
detailed description of the data acquisition and reduc- 
tion may be found there. In brief, OMCl was observed 
using the CFHT with GriF on t he night of Decem ber 
5th 2000. The GriF instrument (|Clenet et alJl2002l) . is 
a combination of the PUEO adaptive optics (AO) sys- 
tem on the CFHT with interferometric spectral scanning. 
The interferometer, a Queensgate ET50WF Fabry-Perot 
(FP), affords a measured spectral resolution A/AA~ 2030, 
that is, 150 kms~^. The detector has a field of view of 
36" X 36" with a pixel scale of 0"035. The region of OMCl 
observed, shown in Fig. ^ consists of four overlapping 
fields and the entire region is centred 15"N and 15"W 
of TCC0016 (05''35™14.''91, -05°22'39"31, J2000). 

Each field has been scanned through the H2 line from 
a wavelength in the far blue wing to a wavelength in the 
far red wing, using step sizes between 4.4 x 10^'*/im and 
4.6 X 10~^/im (^ 65kms~^), allowing adequate sampling 
of the instrumental profile. To prevent superposition of 
different FP orders during scanning, a II2 v=l-0 S(l) in- 
terference filter with a 2.122/^m central wavelength and a 
bandwidth of 0.02/im, was inserted between the FP and 
the detector. For each region and each wavelength, a single 
exposure of 400s was performed. Data reduction includes 
dark and bias subtraction, flat-fielding, bad pixel rejec- 
tion, image recentering, 2D wavelength correction and 
subtraction of sky background obtained from an image 
in the far wing of the II2 profile. Various reference stars 
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Fig. 1. Velocity integrated emission in the S(l) v= 1-0 H2 
emission line covering the full observed field in OMCl. (0,0) 
marks the position of TCC0016 05''35"14!91, -05°22'39'.'31 
(J2000). The star marks the position of BN. Units are arc- 
seconds, where 1 arcsecond = 460 AU. 



for the AO were used: TCC0016 {my ^ 14) in the SE and 
SW field, Parenago 1838 (my ~ 8) in NE and Parenago 
1819 (my - 12) in NW. The FWHM of the point-spread- 
function of stars in the region has been examined yielding 
a spatial resolution of 0'.'15. This near diffraction limited 
spatial resolution was maintained throughout all observa- 
tions. In view of the spatial resolution of 0'.'15, all images 
have been smoothed by a moving average of 3 by 3 pixels 
to improve S /N without significantly degrading the spatial 
resolution. 

Radial velocities were found by choosing a specific po- 
sition on the sky and taking a cut through the cube made 
up of the channel maps. This yielded a set of count rates as 
a function of wavelength constit uting a line profile and a 
Lorentzian is fitted to the profile f'Clcn et et alJbOO^) . The 
position of the peak in emission gives the radial velocity. 
The Lorentzian form was chosen because it represents the 
instrumental profile of the Fabry-Pero t. We note that the 
instrument profile is highly symmetric fjClenet et al .I2OO2I) 
and thus a choice for example of a centroid rather than 
the peak of a Lorentzian is immaterial to an estimation 
of the velocity. We also note that less than 2% of profiles 
are double peaked: the velocity of the brighter of the two 
peaks is chosen when these cases are encountered (Paper 
I)- ^ 

Performing these operations pixel by pixel for the full 
field of view results in very accurate fits for the brighter re- 
gions. Including statistical error in the count rates, relative 
peak wavelengths may be found in bright local zones to 2 
to 3 parts per 10^(3ct), that is, to better than ±lkms^^. 
The quality of the fit decreases with decreasing bright- 
ness. Since an assessment of errors is important in dis- 
cussing the statistical properties of the velocity field, the 
accuracy of determination of the central wavelength has 



been estimated by performing a large number of fits for 
a range of pixel brightness. An empirical correspondence 
between brightness and error in velocity has thereby been 
determined. This may be expressed as follows: 

'2\ 



0.15 + 18.1 exp(-//4. 
2.15exp(-//1.2)) 



8-10" 



(1) 



where cr is the standard deviation of the radial velocity 
in kms^^ and 1 is the brightness in counts • pixel^^ • s^^ 
for the peak emission associated with the Lorentzian fit. 
All pixels with brightness below 0.05 counts • pixel^^ • 
s^^, that is 2% of the maximum brightness, have been 
excluded. Thus the largest uncertainty, a, in the radial 
velocity that may be encountered in our data is 8.6 kms~^ . 

In addition to random errors in the velocity, systematic 
errors may also occur in establishing velocity differences 
between regions which are physically remote. This may 
take place due to distortions in the 2D wavelength correc- 
tion arising from possible mechanical instabilities or drift. 
In this connection, observations of a restricted part of the 
full field comprising two 36"x36"fields to the SE and NW 
were performed in January 2003. Comparisons have been 
made of the velocity fields in the 2000 data, used here, 
and in the 2003 data. Average velocities have been esti- 
mated within regions of about 50x50 pixels (l'.'75x 1'.'75) 
separated by up to 60", the full size of the image. We find 
that velocity differences in the 2000 data and the 2003 
data agree to better than ±lkms~^. Since these results 
were obtained using independent calibration data, they 
provide a clear indication that systematic errors in veloc- 
ity differences between remote regions are not present. 

Conversion of count rates to absolute values of surface 
brightn ess was obtained by comparison with calibrated 
data in lVannier et al.l 1 200 itl . Thus the count rate of 2.15 
counts • pixel"^ • s~^, at the peak of emission to the SE of 
TCC0016, corresponds to a velocity integrated brightness 
of 3.0±0.15 • 10-5 Wm-^sr-i. 

The composite field observed within OMCl, shown in 
Fig. n is 90" X 70". This allows us to address statistical 
issues involving scales ranging from 70 to 3 TO"* AU or 
3.4-10-'^ to 0.15 pc. 

Velocities are expressed as vi^r by assigning a mean 
velocity of a ll the data of 12zfc6 kms"^ in the local stan- 
dard of rest (|Chrvsostomon et aljFmflj ISalas et alJlTaQi 
lO'DellllioOlh . Since we are essentially concerned only with 
the relative velocities of regions within OMCl, the uncer- 
tainty in the value of visr is not material to our discussion. 

3. H2 emission and excitation in the inner region 
of OMCl 

There is a wealth of data for the dense star-forming 
region of OMCl, even if consideration is restricted 
to IR observations of H2 emission. Reference to these 
extensive data, covering features such as the well- 
known "bullets" and "fingers" of H2 emission, and 
very high spatial resolution observations of II2 per- 
formed using the HST, VLT and CFHT may be found 
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data specifically concerned with radial motion measured in 
Orion using H2 IR emission as a tracer dSueai et al.lll995t 
IChrvsostomou et alJ Il997t ISalas et all Il999(l . Excitation 
mechanisms of H2 are however briefly considered here 
since these have a bearing on later discussion. 

The high brightness of H2 emission in the v=l-0 S(l) 
line, frequently exceeding 10~^Wm~^sr^^, indicates that 
the gas is of high column density. Since brightly emit- 
ting clumps are typically 1" to 2" in angular size (1" = 
2.2-10~^pc), the inference is that number densities are 
also high. The most feasible general mechanism of H2 
excitation, for high brightness, is through C-type (mag- 
netic ) shoc ks l|Smith fc Brandl |1990': 'Kauf man fc Neufeldl 
199fia.b.: iTimmermannI ll99Sl"lwilpenbus_et a,l.l l2000t_ 



lHHha.hi: I nmmermanni IIHH^ iWilgenhus et al.l l^niMl 
Vannier et al.ll200lt iLe Boyrlot et alJl2002l: iFlower et all 



2003HKristensen et al .'2003'. '2005^. Low velocity shock ex 



citation in dense gas in regions of OMCl s hown in Fig.U 
has be en anal ysed in detail for e xampl e in IVannier et all 
lj200lh and in iKristensen et al.l l)200.^ . Densities in the 
post-shock gas exceed lO^'cm"''. In addition, the data 
of Paper I graphically illustrate the presence of shocks 
through the clear spatial association of bright H2 emis- 
sion with gas flows. 

In addition to shocks, photon induced excitation of H2, 
in so-called photon dominated regions (PDRs), can play 
a role in yielding the obser ved H2 IR emiss ion. This has 
been discussed in detail in IKristensen et al . (2003) who 
show that PDRs may contribute significantly in specific 
regions, for example in the SE (Peak 2) region, especially 
at the edge of clumps of material in this region. In strongly 
emitting zones, PDRs cannot vield sufficient b rightness, 
even including effects of advection l|Lemaire et al...l999^) . 
high density, and the high UV flux fr om 8^0ri-C, in excess 
of 10 ^ times the standard UV field ijStorzer fc HollenbachI 

4. Statistical Analysis 

4.1. Filtering of data 

The presence of large scale velocity differences in OMC l 
has been known fo r man v ve ars. ISugai et all l|l99,'j) . 
IChrvsostomou et"al] l)l99 7^ and 'Salas et al.l l)l999') have 
all used FP interferometry in OMCl in the IR and stud- 
ied the large scale velocity structure of the hot H2 com- 
ponent. From our observations (Paper I) we find that the 
gas in Peak 1 (NW of BN) is on average 10 kms~^ more 
blueshifted than in Peak 2 (SE of BN: Fig. [H see also 
Fig. 2 in Paper I). This velocity difference has generally 
been viewed in the context of an outflow mechanism from 
the IRc2-comple x and has previously been explained by an 
expanding sheU ('Scovillc ct al."l982': 'Sugai et alJll99 5'). a 



cal wind ijSalas et al.lll999() . The velocity difference could 
however be a manifestation of turbulent energy injected 
into the system and stirring of the gas in large scale vor- 
tices. Thus the nature of the velocity difference remains 
uncertain and in the following we will explore this question 

further. 

\ ■ discussed in iMiesch fc Ballvl (Il994l) . 
iMiesch fc Scalol l|l995|) and iMiesch et alJ l)l999(l . large 
scale trends - such as a general velocity gradient - could 
dominate any statistical quantifier and give misleading 
results. This is due to the fact that the statistical 
methods lose spatial information and therefore depend 
on an absence of systematic trends within the data. 
It may therefore be argued that any large scale trends 
should be removed and only local v elocity fluctuations 
studied. lOssenkopf fc Mac Lowl l(2fl02l) however point out 
that large scale systematic motions may be part of the 
turbulent cascade if these inject energy into the system 
and should therefore only be removed if the turbulence 
studied is exclusively driven on smaller scales. 

In this connection, observations of molecular clouds 
seem to indicate that interstellar turbulence is driven on 
scales larger than the clouds themselves, and arguments 
have been presented that supernovae explosions domi- 
nate (|Mac Low fc Klessenir2004) . However, contributions 
to the continuous driving o f the turbulence ma^ also arise 
from protostellar outflows (jMatzner fc McK'ee 2000") and 
stella r winds from massive O stars (Ma c Low fc KlessenI 
I2OO4I) . OMCl is itself a region of highly active star for- 



mation. Observations ranging from X-ray ( Garmire et alJ 
2000t iFeigelson et al.l [2002.) to radio wave lengths (e.g. 
Churchwell et all Il987t IZanata et all 12004) reveal the 



presence of protostar s, outflows (e .g. 'Gcnzel fc Stutzki' 
il989.) . HH objects Ce.g. lO'Dell et al.ll997.: Doi et al..,2004.) 
and rnainlin e OH, H^^O and SiO r naser emission (e.g . 
'Norris' '1984; 'Mnnton fc Il99,4 iGanime et all h 99j4 

Greenhill et al. . 200 4b). noting that maser emission of the 
nature observed in OMCl is associated with the presence 
of 0-stars. Thus, in OMCl the turbulence could be driven 
locally by the powerful outflow from the BN-IRc2 com- 



.. ^^^^^^ ^ 

as source 1, s 


ource n and BN itself ( 


Menten & Reid"l995'; 


Gezari et al...l998£ iDoeleman et al.l 


1999; 


.Greenhill et alj 


2004al IShuDing et all l2004|). Turbulence mav also be 



bipolar outflow 



or a spheri- 



driven on smaller scales by low mass star formation. 

The alternative mechanism remains, as mentioned ear- 
lier, that the larger scale dynamics is predominantly part 
of a turbulent cascade driven on greater scales. For this 
reason we test the statistical quantifiers explored in this 
paper with the full velocity field and also the residual ve- 
locity fluctuations resulting from removing the large scale 
velocity trends. 

The proce dure of removal o f large scale gradients, or fil- 
tering, follows IMiesch fc Ballvl l)l99-d) . The removal is per- 
formed by convolving the data with a smoothing function 
and then subtracting the smoothed image from the origi- 
nal to obtain the residual image of the small scale velocity 
fluctuations. The residual image should then reveal any 
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small scale structures that were superposed on the large 
scale trends in the original. The filtering or smoothing 
function used here is a two-dimensional equally weighted 
moving average in the form of a square box. The optimal 
size of the filter is the broadest possible that removes the 
large scale gradient but leaves all other features. The 2D 
autocorrelation function of a velocity map can be used for 
detecting velocity gradients. These show up as anticorrela- 
tions in the direction of the g radient and recorrelations i n 
the direction at right angles ijSpicker fc Feitzingerlll988() . 

In the autocorrelation function of our velocity map the 
correlation persists in the NE-SW direction and decorrela- 
tion occur in the NW-SE direction. This indicates a large 
scale velocity gradient in the NW-SE direction which is 
consistent with the velocity difference between Peak 1 and 
2. The widest filter where no anticorrelated sidelobes are 
discernible in the autocorrelation function of the residual 
image is the optimal. By calculating the autocorrelation 
function for filters of varying size, the optimal filter size 
was found to be 14" x 14" (6400 x 6400 AU). All further 
analysis has been carried out on both the original full ve- 
locity field and on the residual image obtained from fil- 
tering with the optimal filter. This latter will be referred 
to as the filtered velocity image or field. We find that it 
is certainly unwise to remove apparent large scale velocity 
trends and that the filtered image is not representative of 
the dynamics in OMCl, as discussed immediately below 
in Sect. 1321 

4.2. Size-linewidth relation 

iLarsonl l)l98l|) identified an empirical relation between 
linewidth (or velocity dispersion) and size for molecular 
clouds and clumps within clouds over scales of ~ 0.1-100 
pc. Larson obtained a power law 



Avobs oc i?" a ~ 0.38 



(2) 



where Avobs is the observed linewidth and R is the effec- 
tive radius of the object. This relationship essentially sup- 
ports a model in which larger motions are associated with 
larger scales. More recent studies give a range of values 
for ct over dimensions of 0.02-100 pc fe.g. iCaselli fc Mverj 
llflfl.'iUPeng et al .119981) showing a varying between 0.2 and 
0.7. Some studies indicate that massive star formi ng re- 
gions tend toward the lower values (e.g. ICaselli fc Mverj 
I1995I) . There are however indications that the relation in 
eq. (2 ) breaks do wn in the most massive star forming re- 
gions JPlume et al. 1997) . A co mprehensive overview was 



given in iGoodman et ai.l l)1998|) . Size-linewidth measure- 
ments require that objects are well defined within a map 
and that a linewidth and size can be objectively assigned 
to objects within the field of observation. The definition 
of objects is obvious for isolated clouds, but encounters a 
certain degree of arbitrariness when dealing with clumps 
in a cloud. 

The structure of OMCl, measured in vibrationally ex- 
cited H2, reveals a very clumpy environment where the 



regions of bright emission are in many circumstances im- 
posed on extended weaker emission. OMCl is thus a case 
study in the difficulty of determining the physical ex- 
tent of individual clumps in a turbulent environment. In 
order to circumvent this problem in computing a size- 
linewidth relation we f ollow the method developed by 
lOssenkopf fc Mac Lowl l)2002() . We remind the reader at 
this stage that our observations of linewidth are indirect. 
That is, we obtain the velocity associated with any pixel 
from the peak of the line profile observed with the low 
resolution Fabry-Perot, see Sect. [3 Thus the linewidth 
associated with any chosen assembly of pixels is given by 
the velocity dispersion of these peak velocities. This is 
in contrast to the standard technique in radio observa- 
tions in which the velocity resolution of observations is a 
fraction of the linewidth, and the linewidth may then be 
obtained directly from the measured lineshape. The tech- 
nique which we use here, that of taking the peak of the 
line fo r each pixel, was also used bv lOssenkopf fc Mac Lowl 
l(2002l) (using centroid velocities) . For this reason we refer 
below to the Larson relationship as a size- velocity disper- 
sion relation, rather than the more familiar size-linewidth. 
In the present case we note that H2 emission is optically 
thin and there are n o optical depth effects, such as those 
discussed in detail in lOssenkopf fc Mac Lowl 1 20021) . 

The average brightness weighted velocity dispersion, 
Avobs, within regions of varying size, is estim ated follow- 
ing the prescription of lOssenkopf fc Mac Lowi (,2002.) . We 
compute the brightness weighted probability distribution 
function of the peak velocities in a circular region with a 
given radius and find the corresponding velocity disper- 
sion from a fit to a gaussian. This is repeated throughout 
the map and for a number of radii, R (see Eq. ||2J)). For 
each radius the average velocity dispersion is calculated 
using intensity weighting. The relation between velocity 
dispersion and size is shown in Fig. El for both the full 
velocity image (upper frame) and the filtered velocity im- 
age (centre frame). Errors shown are the statistical stan- 
dard error on the mean. Within the errors the relation 
for the full velocity field is consistent with a single power 
law over more than two orders of magnitude in R with 
an exponent a — 0.205 ± 0.002, although there also ap- 
pears to be some deviation above 8000 AU. The value of 
the exponent agrees with th e average value of 0.21 ± 0.03 
that lCaselli fc Mverj ^ll995^ found for massive cores in the 
Orion A and B clouds at scales 0.03 - Ipc. It is striking 
that the Larson relationship appears both to hold in this 
very different regime and to have a very similar exponent 
as for massive cores, although we are dealing here with 
highly excited material. 

The appearance of the region south-west of BN with 
fewer and more isolated clumps of emission, suggests 
that the dominating physical processes in this region 
are different from those in Peaks 1 and 2. In Peaks 1 
and 2 the emission is more homogeneous and more spa- 
tially concentrated. All clumps with measura ble radial ve- 
locities south-west of BN are blue-shifted. iNissen et alJ 
((2005.) provide clear evidence that these objects form 
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Fig. 2. Larson size-linewidth relations, Eq. Q. Top: Velocity 
dispersion as a function of size for the full velocity field. A 
power law fit of index 0.205 is overlaid. Centre: Velocity dis- 
persion for the filtered velocity image. Power laws of index 
0.210 and -0.004 are overlaid. Bottom: Velocity dispersion for 
the outflow region (upper line) and Peak 1 [lower line) with 
power law fits overlaid. 



the IR counterpart of an outflow detected in radio ob- 
servations of SiO mascrs associated with a buried O- 
star within OMCl jMc nten fc Rcid 1995; Docl eman ct al.1 
Il999t ICxreenhill et 3111200431^ : IShuping et aLll2004D . For 
this reason we refer to this zone as the "outflow region". 

We have computed the size- velocity dispersion relation 
for the outflow region and Peaks 1 and 2 separately, in each 
case for the full velocity field (that is, without filtering 
of the data - see Sect. I4.1|l . The resulting size- velocity 
dispersion relationships, shown for the outflow region and 
Peak 1, with Peak 2 omitted for clarity, are shown in the 
third frame of Fig. |21 This illustrates that the structure 
adheres to the Larson relationship in the outflow region 
south-west of BN but with a significantly lower slope and a 
somewhat higher proportionality parameter k in Avobs = 
kR" than for Peaks 1 and 2. In the latter two regions k's 
and the exponents a are the same within observational 
error. The lower value of a asso ciated with the outflow 
region supports the conclusion of ICaselli fc Mver^ l)l995|) 
that lower exponents are characteristic of massive star- 
forming regions. 

The relation based on the filtered velocity image in the 
centre frame of Fig. El shows that the velocity dispersion 
is not well represented by a single power law. At radii 



smaller than ~ 1600 AU the relation follows a power law 
of index 0.210, essentially the same as for the full veloc- 
ity field. This similarity of behaviour is expected since the 
filter applied to remove the velocity gradient has a width 
of 14"-- 6400 AU (Sect. Ol), and an equivalent radius 
of 3200AU, and thus should not affect smaller scales. In 
contrast to the full velocity field, the velocity dispersion 
for the filtered data at larger scales is constant within 
observational error as a function of size, in marked con- 
tradiction to the accepted form of the Larson relationship, 
with an index at scales larger than 1600 AU of -0.004 ± 
0.009. In a turbulent medium, where the velocity distribu- 
tion is expected to get broader with increasing size of the 
region, this is unnatural and can only be caused if some of 
the turbulent velocities have been artificially removed. We 
conclude that by removing the large scale gradient we have 
removed some of the turbulent energy. It follows therefore 
that the large scale motions in OMCl , reported both here 
and by a number of other authors JScoville et alJll982l: 
Sugai et 3,1.1 ri99i IChrvsostomou et a1.lll997t ISalas et all 
19991) should be seen as representing real scales of a tur- 
bulent cascade. By implication the turbulence is driven at 
large scales and this suggests that turbulence on the scale 
of the map could be either the result of an energy cascade 
from still greater scales or injection of turbulent energy 
by the large scale outflows associated with massive star 
formation in the IRc2 region. 

This in turn implies that the turbulence in the OMCl 
region does not have a strong injection of energy at scales 
of less than 0. 1 pc and is thus not primarily driven on small 
scales, for example by low mass protostellar outflows. In 
support of this, the energy injected by a high mass stellar 
outflow considerably exceeds that of the sui n total of low 
mass outflows. Elsewhere we have estimated ijNissen et alJ 
^005) that the mass outflow rate is of the order of > 10"'^ 
M0/yr in the blue-shifted outflow, SW of BN, with an 
average velocity of 18 kms~^. Low mass protostellar out- 
flow rates are typic ally three orders of magnitude lower 
l|Richer et al.ll2000|) . with a similar velocity. This simple 
argument on the basis of energetics supports our conclu- 
sion that the injection of energy at the 0.1 pc scale of mas- 
sive stars outweighs on the global scale the energy input 
from low mass stars. However, as we find in Sect.O locally 
as opposed to globally, the character of the turbulence is 
strongly affected by low mass protostellar outflows. 

The velocity gradient, as identifled above, should con- 
sequently not be removed prior to analysis. In the follow- 
ing we thus mainly focus on the full velocity held contain- 
ing velocities on all scales, but for completeness we have 
also included analysis of the filtered velocity image. 



4.3. Peak Velocity PDFs 

The probability distribution function (PDF) of veloci- 
ties is defined in terms of 3D velocity components as 
the number of pixels at a certain velocity vs. veloc- 
ity. Here we are restricted to a set of radial velocities 
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and we estimate the PDF of these velocities sampled 
over the spatial region in the plane of the sky shown 
in Fig. This c orresponds to the PDF of centroid 
velocities used in iMiesch fc Scald (119951) : | Miesch et al.l 
l)l999(l : IOssenkoDf fc Mac Lowll|200j) . An alternative pro- 
cedure is to use the line profile of an optically thin 
line a s an estimate of the PDF al o ng the line-of- 
sight |Falgarone fc PhillipsI Il990l Il99lt iFakarone et^ 
\im4 lOssenkonf & Mac Lowl l2002l) . taking advantage 
of the high velocity resolution of radio-observations. 
lOssenkopf fc Mac Lowl (|2002) provides a comparison of 
the two methods. 

The shape of the wings of the PDF is diagnostic of in- 
termittency, in which, at random scales and/or time inter- 
vals, energy is dissipated as heat in the turbulent medium. 
At any scale the removal of energy must involve those 
elements of the medium with the most prevalent veloci- 
ties. Thus it is the velocities at the centre of the distri- 
bution which becomes most depleted. Thus the velocity 
distribution loses its initial gaussian form, with the centre 
depressed and the wings relatively enhanced. Increasing 
degrees of intermittency create a transition from gaus- 
sian (a exp(— x^)) to exponential (oc exp(— x)) wings in a 
PDF. 

Gaussian PDFs may b e found in studies of decay 



ing supersonic turbulence ijOssenkopf 



fc Mac Lowl 
iBatchelon 



2M 



and incompressible tu rbulence (e.g. iBatcheloJ 11956 
I.Tavesh fc Warhaftlll99ll) . Exponential PDFs can be found 
in the literatur e of simulations, e.g. fro m inelastic colli- 
sions of c louds teicotti fc Fcrrara"2002'l and interactions 
of shells IjChappell fc Scalo 2001.'1 . Further relevant theo- 
retical studies involve the stretched ex ponential form (ex 
exp(— a|x|^)) where is fractional (e.g. lFrisch fc Sornettd 
ll997HEggers fc Wana.1998} . where f3 = 2 reproduces the 
gaussian PDF and /? = 1 the exponential form. Fractional 
P values can arise from random processes whose overall ef- 
fect accumulate in a multiplicative manner (see Sect. l4.4jl . 
Other studies are concerned with wings of PDFs in the 
form of power laws for example arising from stellar winds 
or outflows (,Silk.l995|) . 

U sing line profile dat a for CO f r om various iso- 
topes iFalgarone fc Phillips' ('l990', 'l99lj), iFakarone et alJ 
l)l994 : .Ossenkopf fc Mac Low (,2003) do not find 
simple gaussian behaviour for m ost observations. 
However. I Falgarone fc PhillipsI l)l99l|) showed that most 
of the PDFs could be represented by two gaussians, 
with the wing component being about 3 times broader 
than the core component. The PDFs in the present 
work cannot be so represented. Other work uses centroid 
velocities, co rresponding to the pea k velocities used here . 
For example, iMiesch fc Scalol l|l995l) . lMiesch et all l|l999t) 
found exponential tail s in m ost of their PDFs, while 
lOssenkopf fc Mac Lowl l|2n02l) found that the PDFs of 
the Polaris Flare could be repro duced by two gaussians 
as in lFalgarone~fc PhiUipj l)l99l|) . 

In the present work, the normalized PDF of peak ve- 
locities has been calculated by binning the velocities in 
intervals of lkms~^. The PDF can be estimated by as- 
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Fig. 3. Probability distribution functions of velocities. Top: 
(-I-) PDF for the full velocity field. (O) PDF with removal of 
velocities associated with the clump in Fig. |1| Bottom: PDF 
for the filtered data (see text). Stretched exponentials with 
/3 = 1.10 and /3 = 0.86 respectively are shown for comparison. 
The dashed lines are gaussian fits to the core. 



signing the same weight to every pixel or by weighting 
every velocity with the corresponding brightness in that 
pixel. By testing both methods we found that the chosen 
method of weighting does not influence the shape of the 
PDF significantly. Brightness weighting is less influenced 
by observational noise and is therefore used in the follow- 
ing analysis. 

Fig. shows the peak velocity PDF of the full veloc- 
ity field and of the filtered velocity image. The PDF of 
the filtered velocity image is artificially centered around 
kms~^ through the action of filtering. The PDFs are 
shown on a log-linear scale where a gaussian distribution 
would form a parabola and an exponential a straight line. 
To test the influence of uncertainties in the determina- 
tion of velocities on the shape of the PDF, we have added 
to each velocity in the original dataset a random velocity 
from a gaussian distribution with standard deviation cor- 
responding to the uncertainty in the velocity of the pixel 
in question, as given by Eq. We created a series of 
such velocity fields and calculated the PDF from each of 
them. The PDFs of velocity fields so generated are effec- 
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0.25 0,45 0,65 0,85 




Table 1. Moments and stretching exponents of the peak ve- 
locity PDF for the full velocity field and the filtered velocity 
field 



Fig. 4. The brightly emitting clump that causes the hump in 
the red wing of the PDF (Fig, |3 + symbols) , The xy plane is 
the plane of the sky, the vertical axis shows the radial velocity 
vtsr- The plane at (x,y,43) shows the boundary of the region of 
pixels ignored in the PDF shown as (O) in Fig, |3 The centre 
of this clump is 6','5 E, l','5 N of TCC0016, The colour scale 
refers to brightness in counts pixel" ^s~^, see Sect,|5| 



tively indistinguishable from the PDF of the real data set, 
reflecting the fact that the large number of pixels used 
in the PDF, ~ 3,2-10®, reduces the effect of the observa- 
tional noise to a negligible level. Again, since every bin of 
the histogram is populated by a large number of pixels, the 
relative error, ^/n/n, becomes very small, including tak- 
ing account of the effective 3x3 pixels resolution. Except 
for the very least populated bins the statistical error bars 
would be smaller than the symbols in Fig, |21 and have 
therefore been omitted. 

The PDFs in Fig, O are clearly not well represented 
by gaussians except in the inner core, where gaussians are 
shown as dashed lines in Fig, |31 Our data - for hot dense 
gas - show a dif ferent result to the data rep orted for much 
larger scales in iFal earone fc Phillips' ('1991') (for cool, dif- 
fuse gas and obtained using line profiles), where, as noted, 
a combination of two gaussians fit the observations. The 
wings in the present work can be fitted by stretched expo- 
nentials (full lines). These lines are derived for the full ve- 
locity field data by fitting between +12 and -32 kms~^and 
yield /3=1,10±0,06, A similar fit for the filtered data yields 
/?=0,86±0,05, 



Field 


std-dev / kms ^ 


skewness 


kurtosis 




full 


10.3 


-0,20 


4,4 


1,10 


filtered 


7,8 


-0,042 


6,4 


0,86 



It is clear from Fig, that the PDFs are not sym- 
metric and that the red wing displays large deviations 
from a smooth behaviour, especially at velocities of 50-60 
kms~^, where a hump is seen in the PDF, By inspection 
of the velocity data we found that the hump results from 
a single structure in the field. This structure is the only 
clump of emission with velocities consistently larger than 
50 kms^^ and is shown in Fig, 01 in xy-velocity space. If 
we ignore all pixels in this structure with velocities larger 
than 43 kms~^, then the PDF for the full field no longer 
contains a secondary hump at 55 kms""'^ (diamonds in the 
upper frame in Fig, O)), This suggests that the structure 
involved is an independent entity which does not fit into 
the overall turbulent cascade. The nature of this dense 
fast moving object remains mysterious. 

The shape of the PDF(s) in Fig,|nican be further quan- 
tified by the brightness weighted statistical moments, com- 
puted directly from the dataset as: 



I{x,y)v{x,y) 



std.dev. 
skewness 
kurtosis 



EmapHx,y)[vix,y) ~ 
V ^rnap Hx, y) 

J2mapHx,y)[vix,y) - 

'^^Y.rnapli^-.y) 
'EmapHx,y)[vix,y) - fi]^ 



'^'^EmapHx,y) 



(3) 

(4) 
(5) 
(6) 



where the summations represented by map encompass all 
pixels in the data set. 

The standard deviation quantifies the spread of the 
PDF, the skewness is a measure of the asymmetry and the 
kurtosis characterizes the deviation from a gaussian pro- 
file, A gaussian distribution has a kurtosis of 3 and larger 
values imply that the PDF has relatively more prominent 
wings. An exponential distribution has a kurtosis of 6, 
The calculated values for the PDF from our present data 
are listed in Tabled showing a departure from gaussian 
and from pure exponential behaviour. We conclude that 
the data show clear evidence of intcrmittency. There is 
also a skewness towards the blue. This may arise through 
preferential obscuration of red-shifted flows, compared to 
blue-shifted, in this dusty and partly obscured region. 

4.4. Variance and kurtosis of velocity differences 

The PDFs presented above carry no spatial information. 
In order to retain some of the spatial characteristics of the 
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velocity field and to quantify liow velocities are spatially 
related within the medium, we now construct the proba- 
bility distribution function of velocity differences between 
points separated by a certain distance in the plane of the 
sky, the lag. The resulting distributions should provide a 
more searching test of theoretical models than PDFs of 
peak or centroid velocities. 

Velocity differences, Av — v{r) — v{r — t), are 
used in this analysis. Here both r and r, the lag, are 
two-dimensional space vectors. Previous authors, using 
both observational data and theor etical models, the lat- 
ter for incompressible turbulence ijMiesch fc Scald flOQSt 
iLis et al Miesch et al have described how the 

shape of the PDF changes with different lags. Specifically, 
PDFs were found to exhibit strong non-gaussian forms at 
small lags. 

For homogeneous, isotropic turbulence - as assumed 
here - the same information can be obtained in a more 
compact form by investigating the structure functions as 
a function of lag magnitude, L=|'r|. Structure functions of 
order p are defined as 

Sp{L) =<| v{r) - v{r - r) |f >=<| Av |p> (7) 

This is in fact the traditional manner in which to char- 
acterise turbulence , as Kolmogorov originally prescribed 
( Kolmogorovl Il94ll) . In the inertial range of the tur- 
bulence the structure functions for Kolmogorov turbu- 
lence (non-dissipative, incompressible) obey power laws, 
Sp(L) ~ L''''^^'^ and the exponents are called scaling expo- 
nents. 

Here we study the second order structure function, 
S2(L), that is the variance, and the kurtosis, S4(L)/S2(L)^ 
as a function of lag magnitude, L. In this manner informa- 
tion is obtained concerning the evolution of the PDF as a 
function of lag without the necessity for studying velocity 
difference PDFs directly. We have included a weighting 
function for the velocity differences, see Eqs. (jHl) and El 
and to avoid confusion with the traditional definition of 
the structure functions (Eq. 0) we refer to Eq. ^ as the 
variance function and to Eq. ^ as the kurtosis function. 

Both the variance and the kurtosis functions are pow- 
erful tools when comparing observations with different 
models of turbulence. For example, Kolmogorov turbu- 
lence predicts a power law variation of the variance func- 
tion (see Eq. ijSJ)) with lag ma gnitude, with a scaling expo- 
nent of 2/3, whereas Miesch et al.l l|l999tl find values be- 
tween 0.33 and 1.05 (see below). Numerical and analytical 
studies of driven supersonic magnetohydrodynamic turbu- 
lence &]^_a^_second_order structure function exponent of 
0.74 ((Boldvrev et al.ll2002t) . 

Constructing the variance and kurtosis functions in- 
volves combining data from all pairs of pixels to obtain 
PDFs as a function of lag magnitude, L=|t|. The vari- 
ance and kurtosis of the PDF of velocity differences are 
then 



Ernap T,\t\=l w{r)w{r ~ T){v{r) - v{r - T)f 

(8) 



, . ^ Emap E|x|=L w{i')'w{r ~ T){v{r) - v{r - T)f 

(9) 

where the summations are performed over all pairs of pix- 
els in the whole map that fulfill the requirement |t|=L. 
We have not used L < 4 pixels (0.15 arcsec), since this 
coincides with the spatial resolution of the data. 

In the above, w is some weighting function. An impor- 
tant issue is whether the velocity data should be bright- 
ness weighted (w(r)=I(r) in counts per second) or equally 
weighted (w=l) in calculating the variance and kurtosis 
of the velocity difference PDFs. In order to address this 
problem, we need further to consider errors in the data. 
We have already removed data at the 2% level as described 
in Sect. [5] We now examine whether this cut-off is suffi- 
ciently low for the examination of variance and kurtosis 
functions. We find that a 2% cut-off is too low if we do not 
use intensity weighting but that it is acceptable if intensity 
weighting is included. 

In order to investigate the effect of errors, we use equal 
weighting and calculate the variance function for a number 
of cut-off values in brightness. All the variance functions 
can be approximated by power laws, but we find that the 
scaling exponent, 7 in cf^{L) oc L'' , increases with increas- 
ing cut-off value before converging to a constant value 
when the cut-off value reaches ~ 9% of Imax- This in- 
dicates that pixels with brightness lower than this value 
contribute significantly to the random noise and that they 
should be removed prior to analysis. This analysis empha- 
sizes that caution should always be exercised when equal 
weighting is used since results may show dependence on 
the choice of cut-off. 

The brightness weighted variance function is expected 
to be less influenced by noise. We find that the brightness 
weighted variance function including all pixels is essen- 
tially the same as the equally weighted variance function 
with a brightness cut-off of 9% of Imax ■ Due to the lesser 
influence from noise we use the brightness weighted vari- 
ance and kurtosis functions in the following. 

Errors in the brightness weighted variance function re- 
sulting from uncertainties in the velocities have been cal- 
culated by using the error propagation law and taking the 
Ict uncertainty on the velocity in each pixel from Eq. (|T|l. 
Due to the large number of pixel pairs that goes into the 
calculation of the variance function the relative error is 
typically 10~^ and therefore negligible. 

The variance function is shown in Fig. [S] for the full 
velocity held in a log- log plot. The best fit of the variance 
function to a power law is found to have 7=0.53±0.01 and 
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7 = 0.53 
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Fig. 5. The variance function for the full velocity field. A 
power law form with exponent 0.53 is overlaid. The inset dis- 
plays the variance function of the filtered velocity image com- 
pared to the variance function of the full field. 
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Fig. 6. The kurtosis function, Eq. Q, for the full velocity field 
(solid line). T heoretical predictions for 3 different Reynolds 
numbers from lEggers fc Wane! (ll99Sh are shown for compari- 



this is also shovifn in Fig.|31 For comparison. iMiesch et alJ 
l)l999|) obtained values of 7 between 0.33 and 1.05 for 
sever al molecular clouds at s cales larger than 1.4- 10"^ AU 
while lOssenkopf fc Mac Lowl |2002,) obtained 7 = 0.94 for 
the Polaris Flare at scales larger than 2000 AU. It is note- 
worthy that our value of 0.53±0.01 is significantly lower 
than the Kolmogorov value of 0.67. We also note that 
the variance function is closely related to the Larson rela- 
tionship such that the exponent of the variance function 
should be of the order of twice the Larson exponent. This 
is approximately satisfied here. 

The variance function in Fig. in fact deviates sig- 
nificantly from a single power law, again underlining the 
non-Kolmogorov nature of the gas dynamics. This is par- 
ticularly apparent through a positive deviation around 800 
AU and a negative deviation around 100-200 AU. The cur- 
vature of the structure function may suggest the presence 
of a multi-fractal medium. At all events, several power 
laws operating in different ranges would appear necessary 
to fit the behaviour of the variance function. 

The physical meaning attached to (L) in Fig. El is 
that it approximatel y represents all energies in eddies be- 
low any scale L ( D avidsonI |2004|) . The deviation from a 
power law in Fig. below ~ 2000 AU therefore represents 
more energy at lower scales and implies that there is an 
excess of material in motion at 2000 AU and below. This 
scale suggests that the excess arises from outflow events 
associated with star formation which result in energy in- 
jection at scales below 2000AU. This featur e is apparent 
from t he data presented in Paper 1 and in 'Ni ssen et al.l 
l|2005l) . Conversely the material appears to suffer less ve- 
locity dispersion at scales below 300AU. Thus material 
associated with these smaller scales appears to have dis- 
sipated some of its turbulent energy at larger scales. We 
speculate that these smaller less turbulent scales may be 
associated with protostellar nebulae or perhaps evaporat- 
ing discs around protostars. In earlier work, purely on the 



basis of spatially res olved imaging l|Vannier et al. 1 120011 
iLacombe et"al1l2004j) . a range of preferred scale sizes ly- 
ing between 700 AU and 1100 AU was identified in Peak 
2. In those cases and in the present case we suggest that 
the breaks in power laws are directly linked to scales as- 
sociated with star formation. 

The inset in Fig.jSjshows the effect of the removal of the 
large scale gradient on the variance function. The variance 
function of the filtered velocity image is identical to the 
variance function of the full field for scales smaller than 
half the size of the filter 3000 AU) above which the 
function flattens. It is clear that the filter has no effect 
on very small scales and that the removal of the large 
scale velocity structure causes the variance function to be 
systematically smaller on larger sc ales. This effect o f the 
filtering has also been noted by .Miesch et alJ l)l999|) and 
is most likely an artefact due to the finite size of the map 
llOssenkopf fc Mac Lowll2002|) . 

Turning now to the kurtosis function, the kurtosis is 
a measure of the degree of correlation in the motion at a 
certain scale relative to overall motions seen in the map. 
A velocity field with no spatial correlations would display 
a gaussian shape. As mentioned above, a gaussian has a 
kurtosis of 3. Thus values of kurtosis exceeding 3 indicate 
the strengt h of velocity correlation at sp ecific lags. From 
simulations lOssenkopf fc Mac Lowl l)2002l) found that kur- 
tosis values above 3 can only be verified if the map has 
a scale considerably larger than the lag. As the lag ap- 
proaches the size of the map, the value of the kurtosis 
tends to approach the gaussian value of 3. This arises be- 
cause the map cannot of course contain larger scales than 
the size of the map itself. Thus a kurtosis value around 
3 is always expected at scales approaching the size of the 
map and provides no useful information on the velocity 
correlations at that scale. 

Lo oking specifically at filamentary structure. iLis et alJ 
l)l998i) found non-gaussian wings of velocity difference 
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PDFs at small scales evolving into gaussi ans at larger 
scales , without quoting values of kurtosis. iMiesch et alJ 
1)1999(1 found kurtosis values between 10 and 30 at the 
smallest lags (1.4T0^ - 1.2-10'' AU), with the required 
gaussian value of 3 at large lags of the order of the 
map size (2-10'^ AU). Simila r behaviour was found by 
lOssenkopf fc Mac Lowl l|2002|) . The kurtosis function de- 
rived from our present data is shown in Fig. as a con- 
tinuous line. The kurtosis values approach 30 at small 
lags (100 AU) and decrease gradually until reaching the 
value of 3 at large lags. The shape of the kur tosis function 
shows some similarity to those presented in iMiesch et al.l 
l|l999|) . However a much more striking comparison may 
be found with resu l ts of a mo del of multifractal proc esses 
in lEggers fc Wand l|l998l) . In lEggers fc Wand l|l998() the 
kurtosis, called flatness in that work, is constant out to 
some lag beyond which the kurtosis decreases with in- 
creasing lag. The sha pe of the kurtosis was found by 
lEggers fc Wand l|l998l) to depend on the Reynolds num- 
ber, where R >1.1T0^ marks a transition from a convex 
shape to a concave shape. The shape remains essentially 
unchanged for higher Reynolds numbers. 

In order to compare with results in lEggers fc Wand 



l|1998j) we make an estimate of the Reynolds number of 
the flows in OMCl. The Reynolds number is given by 
Re = ul/v, where u and / are the typical velocity and 
length scale, respectively, associated with the gas flows. 
The dynamical viscosity, v, can be approximated by ;/ ~ 
Awt/i , where A is the mean- free- path of the pa rticles and Vt h 
is the rms thermal velocity (kT//xm^f )i/2 (lFnschlll99fil) . 
Thus, 



Re 



ul 



— ulr 



kT 



(10) 



wher e the cross sect ion, cr, for H2-H2 collisions is 2.7-10 
cm l|Atkinsl Il99 ?l and n is the number density. The 
shocked regions observed here have a typical temperature 
of 2000-3000K and gas densities are 10^ to lO'^ cm'^. The 
velocity u may be approximated by the velocity disper- 
sion in our observed region and u is therefore given by 
twice the standard deviation of the PDF of peak velocities 
(Table The corresponding I is the size of the region ^ 
3-10^^ AU. Using these values, we obtain a Reynolds num- 
ber of Re ~ 3 - 10^. The significance of this value is only 
th at it is very lar ge and greatly exceeds the critical value 
in lEggers fc Wangi. (,1998) for which convex behaviour of 
the kurtosis function vs lag becomes concave. If we were to 
include magnetic effects, the magnetic Reynolds number 
might b e substantially h i gher. 

The lEggers fc Wand l|l998(l description of turbulent 
processes is ad hoc in the sense that it does not result from 
detailed high resolution three dimensional MHD simula- 
tions. Instead some reasonable assumptions are made con- 
cerning the cascade of turbulent energy, assuming it to be 
a stochastic process, isotropic and homogeneous in space, 
with certain additional properties which are briefly sum- 
marized below. Once one assumes that the turbulence is 
isotropic and homogeneous, the turbulence in the Fourier 



domain is represented by a spectrum which is a function 
of a single variable r, corre sponding to a size of t urbu- 
lent eddies. The approach of lEggers fc WangI l|l998l) is to 
impose a recipe for the cascade of energy from larger to 
smaller scales. The recipe states that energy at a given 
scale r is transferred only to a scale r/2. The amount of 
energy transferred, or more precisely the ratio s of the ve- 
locity amplitudes between scale r and r/2 in equilibrium, 
is given by a probability distribution : 



p{s) ^ p5{s - si) + {\ - p)5{s ~ S2) 
with 



(11) 



p = 0.688, 



si = 0.699, 



S2 = 0.947 



where the probability p and parameters si and S2 shown 
above were obtained by fitting to reproduce laboratory 
experimental data. 

Equation Ullfl states that the probability p is high 
(~ 2/3) for the ratio of velocities to be relatively small 
(s — si), which corresponds to a low efficiency of en- 
ergy cascade from one scale to the next. However there 
remains a probability of 1/3 that there is efficient cas- 
cading, s = S2- Starting from some chosen outer scale, this 
recipe allows for a very quick calculation of the energy at 
any given smaller scale by multiplying the probabilities 
for every cascade step lying between these two scales and 
producing the appropriate PDF semi-analytically. Hence 
the term multiplicative turbulence. The viscous cut-off of 
the turbulence is implemented through stopping the fur- 
ther cascading of energy if the Reynolds number at a given 
scale falls below a set value. The velocity fluctuation spec- 
trum is replaced below that scale by an exponentially de- 
caying tail. 

We show in Fig.|H|a compa rison between ou r obser ved 
kurtosis and that presented in lEggers fc Wand l|l998(l for 
Reynolds nu mbers between 1.4 • 10^ and 2.9 • lO'', the 
highest which lEggers fc Wand l|l998(l show, noting once 
more that the form of the kurtosis function appears to 
be independent of Reynolds number above ~ 10^. The 
lower Reynolds number cases are shown purely for com- 
parison. The similarity of the kurtosis function between 
that for our data and that of .Eggers fc Wa ng (J_998) sug- 
gests that the model of a multifractal, multiplic ative tur- 
bulent medium on which 'Egge rs fc Wan is based 
may capture some of the physics of the turbulence in the 
hot component of the gas in OMCl. 

5. Analysis of individual clumps 

The high spatial resolution of the data allows the analysis 
of individual clumps of shocked gas with the same statisti- 
cal methods as above. First, the dimensions of individual 
clumps must be deflned and the following algorithm is 
used to identify the extent of any clump. The data are 
smoothed by a 9 by 9 pixel boxcar, in order that random 
fluctuations become unimportant. A clump is defined as 
a region which encompasses all pixels surrounding a local 
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brightness maximum where the emission is situated on a 
continuously decreasing slope from the maximum bright- 
ness. Thus when we move in any direction starting from 
the pixel of maximum brightness and encounter a pixel 
with higher brightness than the preceding pixel, we define 
this as the boundary of the clump in this direction. As 
a further constraint we only use pixels with a brightness 
larger than 20% of the local brightness maximum. This 
avoids clumps becoming unrealistically large in the outer 
regions with low brightness and few brightness maxima. 
The 20% restriction has no effect in bright congested re- 
gions such as Peak 1 and 2. In passing, we note that we 
cannot rule out that some clumps so identified are the 
result of chance superpositions of two or more isolated 
clumps in the same line-of-sight. We have ignored this 
possibility. 

With the above definition, we have delineated 170 
clumps. These are i n most part e s sentia lly the same fea- 
tures as analyzed in iNissen et al.l l|2005 lh. The number of 
pixels in each clump ranges from 841 to 18026, which cor- 
responds to sizes of clumps between approximately 500 
and 2200 AU. Due to the small size of the clumps, the 
earlier discussion concerning the removal of the large scale 
gradient is irrelevant here, since filtering does not affect 
the velocity field over such restricted regions. 
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5.1. PDFs of dumps 

The PDFs of the 170 individual clumps have been calcu- 
lated. The sensitivity of the shape of the PDF to uncer- 
tainties in the velocities has been investigated using sim- 
ulations in the same way as in Sect. 14. Jl Errors vary from 
clump to clump but do not influence the general shape of 
the PDF in any of the clumps. Eight representative PDFs 
with error estimates are shown in Fig. 13 Error bars shown 
represent the values spanned by the observed data and the 
simulated data. The PDFs for the clumps take on many 
different shapes. For most clumps the shape is complex 
and appears bi- or multimodal. Bi- or multimodal PDFs, 
such as shown in Figs. b, c and d, can result from 
bipolar outflows from one or multiple protostellar objects. 
Most stars form in binary or multiple systems and t he jets 
from these are found to be episodic or pulsed |Larso ^20031 
and references therein). 39 clumps have clear stretched ex- 
ponential wings. Two such clumps are shown in Figs.|7^,f. 
27 clumps have PDFs that can be fitted by a single gaus- 
sian. Two examples are shown in Figs. Efejb. 

The spatial distribution of clumps with gaussian 
PDFs, stretched exponential PDFs and bi- or multimodal 
PDFs is shown in Fig.|Hl Clumps with gaussian PDFs are 
found in all regions of the observed field, except in the 
central region around BN-IRc2, with essentially the same 
distribution as clumps with non-gaussian PDFs. That is, 
there is no tendency for clumps with gaussian PDFs, for 
example, to group together in small areas. Thus it seems 
that the turbulence dominating the BN-KL nebula can 
simultaneously produce both gaussian and non-gaussian 



Fig. 7. PDFs of eight representative clumps, a), b), c) and d) 
are multi-modal, e) and f) are stretched exponential, g) and h) 
are gaussian. Positions are shown in Figs. IH HTUIITTI 

PDFs in clumps with sizes of ~ lO'^ AU, and that the 
clumps with gaussian and non-gaussian PDFs are inter- 
mingled. However there appears to be some mechanism 
that inhibits the formation of clumps with gaussian PDFs 
in the outflow region around BN-IRc2. 

In connection with the above, gaussian PDFs would 
of course arise if the velocity data had a significant noise 
contribution. However it turns out that this is not the case 
save perhaps in 3 of the 27 clumps mentioned above. For 
six of the 27 clumps with gaussian PDFs, the la uncer- 
tainty on the velocities (Eq. 0J) exceeds the half- width 
(one standard deviation) of the gaussian PDF in 1% or 
more of the total number of pixels. The numbers of pixels 
are 1%, 4%, 5%, 28%, 32% and 54% of the totals. In the 
latter three of these six clumps the velocity distribution 
could be dominated by uncertainties in the velocities and 
these data may therefore be spurious, in the sense that 
the gaussian character may be over-emphasized by noise. 

5.2. PDFs of velocity differences: variance functions 
for clumps 

The variance function of the velocity field has also been 
calculated for each clump, using Eq. ©. Most of the vari- 
ance functions appear to be well approximated by power 
laws with a few exceptions, but the value of the scaling 
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Fig. 8. Position of clumps with bi- or multimodal PDF (blue), stretched exponential PDF (green) and gaussian PDF (red). 
The underlying grey background represents the spatial extent of velocity integrated H2 emission at 2.121 nm with brightness 
greater than ~ 2.4 x lO^^Wm^^sr"^. The position of BN is marked with a star, IRc2 with a circle. a,b,c,d,e,f,g,h mark clumps 
shown in Figs. [7| ^ 



exponent 7 varies between essentially zero, that is, a flat 
distribution, and 1.61. The variance functions are shown in 
Fig.Elfor the eight representative clumps shown in Fig.[7| 

The errors in the variance functions have been calcu- 
lated in the same way as for the full region in Sect 14.41 
The magnitude of the errors depend on the brightness (see 
Eq. (^) and the physical size of the clump, where size is 
the dominating parameter. Thus the smaller clumps have 
the larger errors, but even the smallest clumps have rel- 
ative errors of no more than 10%. Two of the clumps in 
Fig. |51 have relative errors of 10% which is of the same 
order as the size of the symbols in Fig.^ Thus the errors 
have been omitted in Fig. |^ Errors in values of 7 arise 
from lack of precision of the power law fit, rather than 
from random errors in the variance. Typical errors in 7 are 
±0.01. Clumps with gaussian PDFs show 7 = 0.0 - 0.97 
including or excluding the 3 clumps with significant noise 
contributions (see Sect. I5.1|l . Clumps with non-gaussian 
PDFs cover the whole range of 7 = 0.0 - 1.61. 

A high value of 7 represents more ordered structure in 
the velocity field and large velocity gradients within the 
clump. This suggests the presence of ordered shock struc- 
ture. Figure ^] shows the spatial distribution of clumps 
with 7 > 1.0, that is, clumps with the highest degree of 
order in the velocity field. The less ordered clumps with 
7 < 0.67, the Kolmogorov value, are shown in Fig. ^2 



Examination of Figs. ^1 and ^2 shows that the value 
of 7 docs not correlate with brightness, that is, high 7 
can occur equally for strong and weak emission clumps. 
Furthermore, in Peaks 1 and 2 - but not in the central 
region around BN-IRc2, see below - clumps with high 7 
may be found mixed with clumps with low 7 in the same 
spatial region. A theory of star formation including tur- 
bulence, self-gravity, bipolar outflows, etc. should thus be 
able to reproduce such a range of exponents of the vari- 
ance function within a limited physical region. 

The majority of clumps with 7 less than the 
Kolmogorov value of 0.67 resides in Peaks 1 and 2 as seen 
in Figs. 1101 and 1111 The major part of these corresponds 
to clumps with measured velocities which display no clear 
structure and whose rela tive velocities with in any clump 
do not exceed 5 kms~^ jNiss en et al.ll2005|) . accordingly 
giving low values of the exponent. We suggest that some 
of these regions may arise from photodissociation regions 
(PDRs) involving an irregular surface containing a vari- 
ety of lines of sight, illu minated by 6 > ^0ri-C , in a model 
essentially as described in lField et all l)l994|) for the PDR 
NGC2023. This would be consistent with a region lacking 
regular structure. However some low 7 clumps are very 
bright and cannot be reconciled with a PDR excitation 
model. These are likely to be shocks travelling across the 
line of sight. 
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Fig. 9. The variance function for the eight individual clumps 
shown in Fig. Q Power law fits are overlaid and the value of 
the exponent is given. 



Many clumps with high 7 (Fig. 1101) are situated 
in the vicinity of the BN-IRc2 complex in the cen- 
tral region. Moving away from BN-IRc2, their occur- 
rence becomes progressively less. Some of the high 7 
clumps possess morphologies which resemble bow shocks 
(e.g. at (-18", 0") and (-21",-7"), see Fig. ^ and 
are part of the blue-shifted lobe of a bipolar out- 
flow from r adio source I, itself associated with a mas- 
sive 0-sta.r llMenten Reidlll 99,4 IPoeleman et a,1.lll999t_ 
lJl2004albl: l^ 



Greenhi ll et alJl2004albl: IShuping et al j2004HNissen et alJ 

2005) . Other high 7 clumps correspond to protostellar 
candidates identified in Paper I and arise from outflows 
created internally in the clumps e.g. at (-30", 35") and (- 
10", 15"), marked A and B in Fig. [101 

In contrast to regions of high 7, regions of low 7 
are notably absent around BN-IRc2 and in the vicin- 
ity of source I, and tend to congregate where regions 
of high 7 are scarce. Our view of turbulent motion, 
based on model s of tu rbulence, like that for example of 
lEggers fc Wand ( 1998lh , deals only with spatial scales and 
not directions. Thus we tend to consider only isotropic 
turbulence. However in the central region between Peaks 1 
and 2, and also to the south-west of BN, the outflow from 
source I imposes a strong constraint on the radial mo- 
tion of clumps in that region, where as mentioned clumps 
show a strong blue-shift. Motion is evidently not isotropic 



here. We suggest that high values of 7 arise in this region 
because effects of turbulence in the radial direction are ef- 
fectively swamped in the variance function by large blue- 
shifted motions and the turbulence tends to 2D rather 
than isotropic 3D. As we move away from the region of 
the blue-shifted outflow, low 7 clumps are once more en- 
countered (see Fig. Ill() . It is an interesting challenge to 
theory to test if imposition of a strong outflow has the 
effects that we observe here on values of 7, as a check on 
our interpretation. Models should also demonstrate the 
absence of gaussian PDFs. 

We now consider high 7 regions outside the blue- 
shifted outflow zone, that is, those which are well-mixed 
with regions of low 7. The formation and evolution of 
protostellar objects involve periods where a high degree 
of order in the velocity field surrounding the protostar is 
expected. A bipolar outfiow encountering the circumstel- 
lar envelope would for example shock-excite the gas while 
creating an organized velocity field, with accompanying 
high 7. The detection of clumps with high values of 7 
could therefore potentially allow an independent method 
of identifying early protostellar objects in the Class 0/1 
phase. This method may prove of value with the advent of 
very high spatial resolution radio maps with the Atacama 
Large Millimetre Array. 

6. Conclusions 

The fraction of gas studied in the present work is highly 
excited and very dense and quite distinct from gas whose 
statistical properties have been studied in earlier work. 
The latter refers to much larger scales and to much cooler, 
lower density and relatively quiescent gas. Nevertheless we 
have found that the hot dense gas in OMCl nicely follows 
the general trends observed in earlier studies, down to the 
smallest, densest scales investigated here. 

Our major conclusions from the statistical analysis of 
flows of H2 in OMCl may be outlin ed as fol l ows: 

1. The size-linewidth relation of iLarsoiil l)l98lj) is re- 
covered at the small scales investigated here. The scaling 
exponent is a = 0.205 ± 0.002, w hich ag:rees with t he av- 
erage value of 0.21 ± 0.03 that ICaselli fc Mver j l|l99,'J) 
found in Orion at scales 0.03 - Ipc, using CO as an indi- 
cator (and which did not in fact include OMCl). 

2. We find that the outflow, shown elsewhere 
to originate from a young and deeply 



l|Nissen et ahl l2r 
buried 0-star, between Peaks 1 and 2, generates structure 
which also follows the Larson relationship. Thus quite dif- 
ferent environments preserve the Larson relationship. The 
velocity dispersion is however higher for any given size 
of object than elsewhere in the held and the exponent is 
signiflcantly smaller than in the full region. 

3. If we use velocity filtered data, removing the largest 
scales which appear as an apparent gradient of velocity, 
the Larson relationship breaks down. We conclude that 
the large scale motions are an inherent part of the turbu- 
lent cascade and should therefore not be removed from the 
data. Large scale injection of turbulent energy is the domi- 
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Fig. 10. Clumps with the variance function scaling exponent 7 > 1.0. As in Fig.|Hl the underlying grey background represents 
the spatial extent of velocity integrated H2 emission at 2.121 fim with brightness greater than ^ 2.4 x 10~®Wm~^sr~^ . The 
position of BN is marked with a star, IRc2 with a circle. Colours represent brightness in counts per pixel per second in the 
clumps analyzed in sect.|3 A and B mark sites of protostellar candidates. a,b,d,f mark clumps shown in Figs. 0121 




Fig. 11. Clumps with the variance function scaling exponent 7 < 0.67. Grey background as in Fig llOl The position of BN is 
marked with a star, IRc2 with a circle. Colours represent brightness in counts per pixel per second in the clumps analyzed in 
Sect.|K| e,h mark clumps shown in Figs. [7| 1^1 
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nant process and outflows from massive star(s) in the IRc2 
complex may contribute a substantial part of the driving 
of the turbulence at the scale of 0.1 pc. 

4. The probability distribution function (PDF) of ve- 
locities is best fitted by an exponential or a weakly 
stretched exponential and departs strongly from gaussian. 
This suggests that the turbulence in the region studied 
here is characterized by intermittency. 

5. The variance function of the velocity differences is 
not well represented by a single power law. A multifrac- 
tal model is implied. The best fit scaling component is 
0.53±0.01, significantly lower than the Kolmogorov value 
of 0.67, underlining the non-Kolmogorov nature of the tur- 
bulence. 

6. The behaviour of the variance function shows that 
there is a preferred scale size in the medium below 2000 
AU reflecting the presence of protostellar zones of this 
dimension and below. There is also clear evidence of less 
turbulent structure below ~300AU, perhaps representing 
a population of protostellar disks which have expelled part 
of their turbulent energy. 

7. As further evidence of the multifractal nature of 
the medium, the kurtosis function of the velocity differ- 
ences closely res embles that of the multifractal model of 
lEggers fc Wa^ ijlQQS,) for high Reynolds numbers. 

8. Analysis of 170 individual clumps with sizes between 
500 and 2200 AU opens a window on a regime of scales 
that has never before been explored using statistical tech- 
niques. Studies at these scales reveal considerable diversity 
of velocity PDFs, with some gaussian, some showing evi- 
dence of intermittency and many with complex structure, 
reflecting multiple outflows. Variance functions are ap- 
proximately power laws, with exponents varying between 
zero and 1.61. There is no spatial association between high 
and low exponent clumps nor is there any association be- 
tween gaussian PDFs and high or low exponents. 

9. To emphasize the last point, clumps with a variety 
of forms of the PDF and different values of the scaling 
exponent of the variance function are found in the same 
spatial region and in the whole region covered by the ob- 
servations. 

10. An outflow region associated with a deeply embed- 
ded 0-star between Peaks 1 and 2 shows markedly differ- 
ent statistical characteristics from Peaks 1 and 2. Clumps 
in the vicinity of the BN-IRc2 complex typically show high 
scaling exponents in the variance function. This may be 
due to the action of the energetic outflow imposed on in- 
trinsic turbulent motion within individual clumps. At all 
events high and low exponent clumps are spatially segre- 
gated in this zone, with low exponent clumps found only 
at the edges. 

11. The outflow region apart, our results suggest that 
analysis using variance functions could be a useful way 
in which to establish the presence of early star-forming 
regions. 

The above results constitute a challenge for numerical 
simulations of turbulence in star forming regions. To our 
knowledge the resolution of any present numerical simu- 



lation is substantially lower than the resolution of these 
observations, but advances in computer technology will 
soon allow simulations to reach the scales encountered 
here. A self-consistent theory of star formation including 
self-gravity and MHD turbulence should be able to repro- 
duce the features outlined in this paper. Thus the size- 
linewidth relation, the probability distribution function of 
peak velocities and the variance and the kurtosis of veloc- 
ity differences as a function of lag form a dataset which 
models of star-forming regions should aim to reproduce. 
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